4 Incompressible Euler
4.1 Derivation
In the previous section, we derived mass and momentum balance according to
\[ \begin{aligned} \partial_t \rho + \nabla \cdot \left( \rho \mathbf v \right) & = 0 \label{chap1::massbalance} \\ \partial_t \left( \rho \mathbf v \right) + \nabla \cdot \mathbf \Pi & = \rho \mathbf b, \label{chap1::mombalance} \end{aligned} \]
with \(\rho \mathbf v\) denoting the flux density vector with components \(\rho v_i\) and \(\mathbf \Pi\) denoting the second order tensor \(\rho \mathbf v \otimes \mathbf v - \mathbf \sigma\) with components \(\Pi_{ij} = \rho v_i v_j - \sigma_{ij}\).
Observation: We now have four equations for 13 unknown \(\rho, v_i, \sigma_{ij}\). Note, that considering the mass balance alone, we have one equation for four unknowns: \(\rho, v_i\), which is also undetermined. Considering another balance law, here the momentum balance, however, doesn’t help to close the problem. We need a different path!
General strategy: Let’s narrow down material properties of the considered medium for the physical regime that we would like our mathematical process model to cover. Such specialization strategies aim at enriching the mathematical model with additional closure relations. We hence loose generality of the model, yet gain a necessary prerequisite for well posedness, hence solvability.
Let’s assume that the Cauchy stress can be written as \[ \mathbf \sigma = - p \mathbf I + \mathbf \tau, \]
in which \(p\) stands for the pressure, hence stresses due to isotropic volume change and \(\tau\) denotes the deviatoric part of \(\mathbf \sigma\). For now, we simply accept the fact that this decomposition makes sense. Note, that a detailed discussion of the Cauchy stress tensor and its decomposition will be done later.
With that decomposition, we can re-write mass and momentum balance as
\[ \begin{aligned} \tfrac{1}{\rho} \tfrac{D}{Dt} \rho + \nabla \cdot \mathbf v & = 0 \\ \tfrac{D}{Dt} \mathbf v & = - \tfrac{1}{\rho} \nabla p + \tfrac{1}{\rho} \nabla \cdot \mathbf \tau + \mathbf b. \end{aligned} \]
Let’s now consider a medium (i) and physical situation (ii) subject to the following two assumptions:
- \(|\frac{1}{\rho} \frac{D}{Dt} \rho | << |\nabla \cdot \mathbf v |\)
- \(|\nabla \cdot \mathbf \tau | << | \nabla p|\).
The first assumption (i) means that the total change in density of the medium is much smaller than the dominant term in \(\nabla \cdot \mathbf v = \partial_x u + \partial_y v + \partial_z w\), whereas the second assumption means that the deviatoric stresses can be neglected in the momentum flux density vector \(\mathbf \Pi\).
We can use this to simplify mass and momentum balance equations. From Equation 3.7 we get
\[ \begin{aligned} 0 = \underbrace{\frac{1}{\rho} \frac{D}{Dt} \rho}_{\rightarrow 0} + \nabla \cdot \mathbf v = \nabla \cdot \mathbf v \end{aligned} \]
We can conclude that under these assumptions we can expect a divergence free velocity field, hence no sinks and sources in the domain. From Equation 3.8 on the other hand, we get
\[ \begin{aligned} \rho \partial_t \mathbf v + \underbrace{\mathbf v \partial_t \rho + \nabla \cdot ( \rho \mathbf v) \mathbf v}_{\mathbf v \left(\partial_t \rho + \nabla \cdot ( \rho \mathbf v) \right) = 0} + \nabla \mathbf v (\rho \mathbf v) &= - \nabla p + \underbrace{\nabla \cdot \mathbf \tau}_{\rightarrow 0} + \rho \mathbf b \\ \Leftrightarrow \qquad \partial_t \mathbf v + \mathbf v \cdot \nabla \mathbf v &= - \frac{1}{\rho}\nabla p + \mathbf b \end{aligned} \]
As the deviatoric stresses drop out of the momentum balance, we are also losing unknowns.
The above assumptions (i) and (ii) are called the ideal, incompressible fluid approximation, the arising mathematical model correspondingly the ideal, incompressible fluid model. It is also referred to as the incompressible Euler equations and concisely reads:
\[ \begin{aligned} \nabla \cdot \mathbf v & = 0\\ \partial_t \mathbf v + \mathbf v \cdot \nabla \mathbf v &= - \frac{1}{\rho}\nabla p + \mathbf b \label{eq:euler} \end{aligned} \tag{4.1}\]
ALternatively, we can write it in index notation:
\[ \begin{aligned} \partial_i v_i & = 0\\ \partial_t v_i + v_j \partial_j v_i &= - \tfrac{1}{\rho}\partial_i p + b_i \end{aligned} \tag{4.2}\]
Assuming we know the constant density \(\rho\) and the body force is also known, for instance when given as the gravitational acceleration \(\mathbf b = \mathbf g\), then the incompressible Euler equations comprise a set of four equations in the four unknowns pressure \(p\) and velocity \(\mathbf v\). We are able to solve it, given physically meaningful initial and boundary conditions.
4.2 Decomposition of the velocity field
Let’s consider the velocity field around an arbitrary point \(\mathbf x\). Expanding the velocity field in a Taylor series yields
\[ \mathbf v \left( \mathbf x + \mathbf r \right) = \mathbf v \left( \mathbf x \right) + \nabla \mathbf v \cdot \mathbf r + \text{higher order terms} \]
Here, \(\nabla \mathbf v\) is a second order tensor with components \(\partial_j v_i\). Following Section 2.4, the velocity gradient can be decomposed into a symmetric and an anti-symmetric part:
\[ \nabla \mathbf v = \mathbf W + \mathbf D, \]
with
\[ \mathbf W = \tfrac{1}{2} \left( \nabla \mathbf v - \nabla \mathbf v^T \right) \quad \text{and} \quad \mathbf D = \tfrac{1}{2} \left( \nabla \mathbf v + \nabla \mathbf v^T \right), \]
or in components
\[ W_{ij} = \tfrac{1}{2} \left( \partial_j v_i - \partial_i v_j \right) \quad \text{and} \quad D_{ij} = \tfrac{1}{2} \left( \partial_j v_i + \partial_i v_j \right). \]
\(\mathbf D\) is called the strain-rate tensor that measures the strain or distorsion of a medium. Note, that the strain-rate tensor is symmetric and can therefore be rotated into a coordinate system, in which it is represented as a diagonal matrix. The diagonal entries \(\partial_i v_i\) account for dilation or compression in the three principal directions.
\(\mathbf W\) on the other hand denotes the spin tensor and accounts for local rotational motion or spin. \(\mathbf W\) is skew-symmetric and has zero entries on the diagonal, hence a vanishing trace. Following Section 2.4, we get
\[ \mathbf W \cdot r = \mathbf w \times \mathbf r \]
with \(\mathbf w\) being the axial vector of \(\mathbf W\) defined as
\[ \mathbf w = \tfrac{1}{2} \text{curl} \, \mathbf v = \tfrac{1}{2} \nabla \times \mathbf v. \]
Recall, that \(\mathbf w \times \mathbf r\) indicates a rigid body rotation of a point \(\mathbf r\) around rotation axis \(\mathbf w = w \mathbf n\). All in all, we get
\[ \underbrace{\mathbf v \left( \mathbf x + \mathbf r \right)}_{\text{velocity field}} = \underbrace{\mathbf v \left( \mathbf x \right)}_{\text{translational motion}} + \underbrace{\mathbf w \times \mathbf r }_{\text{rotational motion}} + \underbrace{\mathbf D \cdot \mathbf r}_{\text{internal deformation}}. \]
4.3 Vorticity dynamics
We now define the so-called vorticity \(\mathbf \omega = \nabla \times \mathbf v\) (which is the axial vector for \(2 \mathbf W\)). Our next goal is to understand its evolution under the ideal, incompressible fluid assumption!
By making use of the vector operator identity
\[ - \left( \mathbf v \cdot \nabla \right) \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) = \mathbf v \times \left( \nabla \times \mathbf v \right) \]
we can write the momentum equation of the Euler equations Equation 4.2 in the so-called Lamb representation (1895):
\[ \partial_t \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right) = \mathbf b - \tfrac{1}{\rho} \nabla p. \]
Applying the curl operation \(\nabla \times\) to both sides of the equation yields
\[ \nabla \times \left( \partial_t \mathbf v \right) + \nabla \times \left(\nabla \left( \tfrac{\mathbf v^2}{2} \right)\right) - \nabla \times \left(\mathbf v \times \left( \nabla \times \mathbf v \right)\right) = \nabla \times \mathbf b - \tfrac{1}{\rho} \nabla \times \nabla p, \]
which simplifies as we identify \(\mathbf \omega = \nabla \times \mathbf v\) and consider that \(\nabla \times \nabla *\) vanishes:
\[ \partial_t \mathbf \omega - \nabla \times \left( \mathbf v \times \mathbf \omega \right) = \nabla \times \mathbf b. \]
We exploit a second vector operator identity
\[ \nabla \times \left( \mathbf a \times \mathbf b \right) = \left( \mathbf b \cdot \nabla \right) \mathbf a - \left( \mathbf a \cdot \nabla \right) \mathbf b - \mathbf a \left( \nabla \cdot \mathbf b \right) - \mathbf b \left( \nabla \cdot \mathbf a \right), \]
which for \(\mathbf a = \mathbf v\) and \(\mathbf b = \mathbf \omega\) and considering \(\nabla \cdot \left(\nabla \times * \right) = 0\) as well as continuity (\(\nabla \cdot \mathbf v=0\)) yields
\[ \partial_t \mathbf \omega = \left( \mathbf \omega \cdot \nabla \right) \mathbf v - \left( \mathbf v \cdot \nabla \right) \mathbf \omega + \nabla \times \mathbf b. \]
Alternatively, written in terms of the material derivative:
\[ \tfrac{D}{Dt} \mathbf \omega = \left( \mathbf \omega \cdot \nabla \right) \mathbf v + \nabla \times \mathbf b. \]
This means that total change of the vorticity evolution is due to vorticity production by means of existing vorticity in a non-uniform velocity field, and a second part due to the vortex components of the mass force field \(\mathbf b\). If the latter constitutes a potential, hence \(\mathbf b = \nabla B\), we get
\[ \tfrac{D}{Dt} \mathbf \omega = \left( \mathbf \omega \cdot \nabla \right) \mathbf v. \]
4.4 Irrotational motion and potential flow
If the mass force is a potential and the initial vorticity is zero everywhere, we get
\[ \tfrac{D}{Dt} \mathbf \omega = 0 \]
always! As a consequence the velocity field is characterized by a zero vorticity, hence \(\nabla \times \mathbf v = 0\). Due to mass conservation, it also has a vanishing divergence \(\nabla \cdot \mathbf v=0\). In such a situation, the velocity itself is given as a potential, meaning there is a scalar function \(\phi\), referred to as the velocity potential, such that
\[ \mathbf v = \nabla \phi \]
In combination with the continuity equation, we get
\[ 0 = \nabla \cdot \mathbf v = \nabla \cdot \nabla \phi = \triangle \phi \]
4.5 Bernoulli integral
We saw before that the Euler equations can be put into Lamb representation. Let’s integrate Lamb representation along a streamline of the flow.
\[ \partial_t \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right) = \mathbf b - \tfrac{1}{\rho} \nabla p. \]
Let’s integrate this expression along a streamline \(\mathbf x(t)\) starting from coordinate \(\mathbf x_1\) up to \(\mathbf x_2\), which for simplicity are denoted by 1 and 2:
\[ \int_1^2 \left(\partial_t \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right)\right) \cdot d \mathbf x = \int_1^2 \mathbf b \cdot d \mathbf x - \int_1^2 \tfrac{1}{\rho} \nabla p \cdot d \mathbf x. \]
If the mass force is once again given as a potential \(\mathbf b = - \nabla U\) and by furthermore introducing the pressure increment along the integration path \(d p = \nabla p \cdot d \mathbf x\), we get
\[ \int_1^2 \left(\partial_t \mathbf v \right) \cdot d \mathbf x + \int_1^2 \nabla \left( \tfrac{\mathbf v^2}{2} + U \right)\cdot d \mathbf x + \int_{p_1}^{p_2} \tfrac{1}{\rho} dp = \int_1^2 \left( \mathbf v \times \left( \nabla \times \mathbf v \right)\right) \cdot d \mathbf x \]
Note that \(\nabla \left( \tfrac{\mathbf v^2}{2} + U \right)\) can be written as the total differential \(\tfrac{d}{dt} \left( \tfrac{\mathbf v^2}{2} + U \right) dt\), such that we end up with a general formulation of the Bernoulli equation:
\[ \int_1^2 \left(\partial_t \mathbf v \right) \cdot d \mathbf x + \left[ \tfrac{\mathbf v^2}{2} + U \right]_1^2 + \int_{p_1}^{p_2} \tfrac{dp}{\rho} = \int_1^2 \left( \mathbf v \times \left( \nabla \times \mathbf v \right)\right) \cdot d \mathbf x \]
Let’s step through the individual terms: For stationary processes, the first term vanishes. Assuming a so-called barotropic fluid, in which pressure is a function of density alone (\(p=p(\rho)\)) and vice versa (\(\rho=\rho(p)\)), we can introduce a pressure function \(P(p):= \int_{p_0}^p \tfrac{d \tilde p}{\rho(\tilde p)}\), which allows to write
\[ \int_{p_1}^{p_2} \tfrac{dp}{\rho} = P(p_2) - P(p_1) =: P_2 - P_1 \]
The right hand side term vanishes whenever the integration path \(d \mathbf x\) is aligned with the particle velocity \(\mathbf v\), such that the above statement is always true along streamlines. Note that if we are facing irrotational flow, such as discussed in the previous section, the RHS contribution vanishes always and the Bernoulli equation holds for all locations in space.
Being even more specific, we assume a medium in a stationary regime, for which \(P=\tfrac{p}{\rho}\), e.g. water and mass forces to be given as gravitational acceleration \(U=gz\) and get a better known form of the Bernoulli equations, which reads
\[ \tfrac{p}{\rho} + \tfrac{v^2}{2} + g z = const \]
Often the static pressure \(p\) is combined with the pressure due to gravitational acceleration into the piezometric pressure \(p^* = p + \rho g z\), hence
\[ p^* + \rho\tfrac{v^2}{2} = const \]
4.6 Bernoulli equation in a rotating reference frame
The momentum equations in Equation 4.1 can be written in a reference frame that rotates at \(\mathbf \Omega\). The Lamb representation is then
\[ \partial_t \mathbf v + \nabla \left( \frac{\mathbf{v}^2}{2} - \frac{\left|\mathbf{\Omega}\right|^2\left|\mathbf{x}\right|^2}{2} \right) + 2 \mathbf{\Omega} \times \mathbf{v} - \mathbf{v} \times \nabla \times \mathbf{v} = - \frac{1}{\rho}\nabla p + \mathbf b. \]
Following similar arguments than in section Section 4.5, we can derive the Bernoulli equation, which for stationary (\(\partial_t \mathbf v =0\)), irrotational flow (\(\nabla \times \mathbf{v} = 0\)) and a fluid at rest with respect to the rotating reference frame is
\[ \frac{p}{\rho} - \frac{\left|\mathbf{\Omega}\right|^2\left|\mathbf{x}\right|^2}{2} + gz = \text{const}. \]
We get
\[ \frac{p_2}{\rho} - \frac{\left|\mathbf{\Omega}\right|^2 x_2^2}{2} + gz_2 = \frac{p_1}{\rho}, \]
with \(p_1=p_2=p_\mathrm{ambient}\), we get
\[ z_2 = \frac{\left|\mathbf{\Omega}\right|^2 x_2^2}{2 g}, \]
and so the surface is curved to form a rotational paraboloid.